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Abstract 

An  important  consideration  for  Autonomous  Underwater  Vehicles  in  shallow  water  is  station  keeping. 
Station  keeping  is  the  ability  of  the  vehicle  to  maintain  position  and  orientation  with  regard  to  a  reference  object.  In 
shallow  water  AUV  operations,  where  large  hydrodynamic  forces  are  developed  due  to  waves,  knowledge  of  the  sea 
is  critical  to  allow  for  the  design  of  a  control  system  that  will  enable  the  vehicle  to  accurately  navigate  and  position 
itself  while  in  the  presence  of  non-deterministic  disturbances. 

Recently  it  has  been  shown  that  it  is  possible  to  measure  the  shape  of  the  sea  surface  using  remote  sensors 
such  as  acoustic  probes,  lasers  or  short  wavelength  radar.  These  measurements  have  made  it  possible  to  develop  so 
called  "predictive"  control  strategies  for  many  surface  applications  including  hydrofoil  operations  and  craning 
operations  between  vessels. 

The  ability  to  develop  predictive  control  strategies  for  underwater  vehicles  is  limited  by  the  ability  to 
measure  the  environmental  disturbances  acting  on  the  vehicle.  In  this  paper  we  describe  the  design  of  a  model  based 
predictive  controller  that  employs  sub-surface  sensors  for  disturbance  prediction,  to  reduce  wave  induced  effects  on 
vehicle  station  keeping.  Using  linear  wave  theory  and  recursive  estimation,  an  Auto  Regressive  (AR)  model  of  the 
sea  spectrum  is  developed.  This  dynamic  model  is  then  used  to  develop  a  forward  predictor/estimator  which  is 
embodied  in  the  controller  to  cancel  the  predictable  portion  of  the  non-deterministic  disturbance  on  the  vehicle 
thereby  minimizing  position  error.  Through  simulation  the  ability  of  the  Naval  Postgraduate  School's  "PHOENIX" 
AUV  to  maintain  longitudinal  position  while  subjected  to  actual  wave  disturbance  data  from  coastal  Monterey  Bay 
is  shown. 


Introduction 


For  an  Autonomous  Underwater  Vehicle  to  operate  with  a  high  degree  of  reliability,  disturbances  and  their 
effects  on  the  AUV  must  be  modeled  mathematically  with  an  adequate  degree  of  accuracy.  The  main  source  of  the 
dynamic  disturbances  encountered  by  underwater  vehicles  or  submerged  vessels  are  wave  and  current  induced. 
These  disturbance  forces  arise  from  buoyant  and  inertial  effects  due  to  ocean  wave  kinematics.  The  effects  of 
incident  waves  on  a  submerged  body  can  be  divided  up  in  several  categories.  The  largest,  the  first  order  forces  act  at 
the  incident  wave  frequency.  These  forces  move  the  vehicle,  but  usually  result  in  oscillations  about  a  mean  state. 
Second  order  forces,  which  are  the  result  of  wave  diffraction  and  wave  interaction,  have  several  different  frequency 
components. 

Wave  diffraction  of  a  single  frequency  wave  results  in  a  steady  force  and  a  varying  force  at  twice  the  wave 
frequency.  The  double  frequency  force  may  be  neglected,  as  the  inertia  of  the  underwater  vehicle  effectively  filters 
it.  Interaction  of  waves  at  different  frequencies  also  results  in  forces.  These  forces  consist  of  a  component  acting  at 
the  sum  of  the  wave  frequencies  and  a  component  acting  at  the  difference  of  the  wave  frequencies.  The  force  caused 
by  the  summation  of  the  wave  frequencies,  again,  may  be  neglected,  as  it  is  also  filtered  by  the  vehicle's  inertia.  The 
difference  frequency  component  results  in  a  slowly  varying  force  on  the  AUV.  To  best  assess  the  performance  of  a 
vehicle  in  shallow  water,  it  is  necessary  to  generate  a  dynamic  model  representing  the  wave  induced  water  velocity 
and  acceleration  to  properly  model  the  wave  disturbance  forces. 

Since  there  is  strong  uncertainty  in  ocean  wave  fields,  in  the  past  a  spectral  approach  has  been  best  suited 
for  describing  the  intensity  and  frequency  characteristics  of  ocean  wave  fields.  During  the  late  1980's,  P.  Spanos 
(1986)  presented  research  on  digital  simulation  of  ocean  wave  disturbances.  His  research  expanded  on  the  use  of 
Linear  Prediction  Theory  (Makhoul,  1975)  using  Autoregressive  (AR)  and  Autoregressive  Moving  Average 
(ARMA)  models.  Using  the  former  spectral  method,  considerable  numerical  difficulties  were  encountered  in  getting 
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realistic  simulations,  while  the  latter  methods  experienced  difficulties  accurately  approximating  a  given  spectrum, 
and  an  accurate  approximation  could  only  be  feasible  when  using  a  high  order  AR  scheme. 

In  this  paper,  seaway  dynamics  will  be  modeled  as  a  linear  Auto  Regressive  system  subjected  to  random 
white  noise  input,  where  the  linear  system  transfer  function  is  approximately  mapped  to  the  local  spectrum.  Using  an 
LQE  approach  a  predictor,  developed  from  the  AR  model,  is  used  to  predict  wave  states  at  future  times  based  on 
present  measurements.  Finally,  this  paper  will  develop  a  control  scheme  to  make  use  of  the  disturbance  estimation 
and  show  that  station  keeping,  in  the  longitudinal  (surge)  direction,  relative  to  an  object,  in  the  presence  of  shallow 
water  waves  is  possible. 


Linear  Wave  Theory  with  Applications  to  Shallow  water 

The  simplest  free  surface  wave  formation,  which  nevertheless  has  great  practical  significance,  is  the  plane 
progressive  wave  system.  This  motion  is  two  dimensional,  sinusoidal  in  time  with  angular  frequency  to,  and 
propagates  with  a  phase  velocity  cp  such  that  to  an  observer  moving  with  this  speed  the  wave  appears  to  be 
stationary.  A  cartesian  coordinate  system  is  adopted,  see  Figure  1,  with  z  =  0  the  plane  of  the  undisturbed  free 
surface  (still  water  level)  and  the  z-axis  positive  upwards.  The  vertical  elevation  of  any  point  on  the  free  surface  may 
be  defined  by  a  function  z  =  T With  these  requirements,  the  free  surface  elevation  must  be  of  the  general  form 

r\(x,t)  =  Acos(kx-(tit),  (1) 

where  the  positive  x-axis  is  chosen  to  coincide  with  the  direction  of  wave  propagation.  Here  A  is  the  wave 
amplitude,  and  the  parameter 
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is  the  wavenumber,  the  number  of  waves  per  unit  distance  along  the  x-axis.  Clearly 
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where  the  wavelength  L  is  the  distance  between  successive  points  on  the  wave  with  the  same  phase.  Figure  1 . 

The  solution  of  this  problem  is  expressed  in  terms  of  a  two  dimensional  velocity  potential,  which  must 
satisfy  Faplace's  equation 


V2(|)=0, 


(4) 


and  appropriate  boundary  conditions.  Furthermore,  4*  must  yield  the  wave  elevation  (1)  from 


1  3<b 

B  = - — - 

g  dt 


(5) 


Equation  (5)  is  the  linearized  dynamic  boundary  condition  on  the  free  surface  and  is  an  expression  of  the  fact, 
through  Bernoulli's  equation,  that  the  pressure  on  the  free  surface  must  be  the  same  as  the  ambient  atmospheric 
pressure.  An  appropriate  boundary  condition  on  the  sea  bottom  is 

|U  at  z  =  —  H  ,  (6) 

az 

i.e.,  the  bottom  at  depth  H  is  a  rigid  impermeable  plane.  Finally,  the  free  surface  boundary  condition  is 


2 


Equation  (7)  is  a  combined  dynamic  and  kinematic  surface  boundary  condition.  The  dynamic  condition  was 
discussed  earlier,  while  the  kinematic  condition  simply  states. 


3r|  _  30 
3 1  3  z 


(8) 


that  the  vertical  velocities  of  the  free  surface  and  fluid  particles  are  the  same.  Combining  (5)  and  (8),  (7)  is  obtained, 
ignoring  the  small  departures  of  the  free  surface  q  from  the  horizontal  orientation  z  =  0. 

From  the  requirements  of  the  problem,  it  is  clear  that  the  velocity  potential  0  must  be  sinusoidal  in  the  same 
sense  as  (1);  therefore  a  solution  of  the  form 

0(x,  z,t)  =  sin(/cx-cot)F(z) ,  (9) 


is  attempted.  Substituting  (9)  into  (4),  F(z)  must  satisfy  the  ordinary  differential  equation 


d2F 

dt2 


- k2F=0 , 


throughout  the  domain  of  the  fluid.  The  solution  to  ( 10)  satisfying  the  bottom  boundary  condition  (6)  is 

F  =  Acosh(k(z  +  H)). 


(10) 


(11) 


Substitution  of  (9)  and  (11)  into  the  surface  boundary  condition  (7)  yields  an  important  relationship  between  the 
wavenumber  k  and  the  frequency  to. 


(S)2  =gk  tanhOW) ,  (12) 

which  is  called  the  dispersion  relationship.  The  surface  elevation  q  follows  from  (5) 

r\(x,t)  =  acos(kx-(£>t) ,  (13) 


with  the  amplitude  a  given  by 


a  =^Lcosh(T/7) . 
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Substitution  of  (1 1 )  and  ( 14)  into  the  velocity  potential  function  (9)  yields 


,  ,  ag  cosh(k(z  +  H))  . 

0(x,  z,  t)  = - 77777. —  sin(Cr-CDf) . 


to 


cosh  (kH) 


(14) 


(15) 


An  underlying  assumption  associated  with  potential  flow  theory  and  Laplace's  Equation,  is  that  the  velocity 
field  can  be  expressed  as  the  gradient  of  a  velocity  potential  function.  This  allows  the  expressions  for  the  velocities 
in  the  x  and  z  direction  to  be  given  as 


(16) 


Using  (16)  and  the  linearized  form  of  the  Bernoulli's  equation  (17), 

3d) 

P  =  ~  P^: — P  gz, 
dt 

the  expressions  for  the  fluid  velocity  and  pressure  fields  are: 


gk  cosh(£(z  +  //)) 

it  =  a - 

co  cosh  (kH) 

gk  sinh (k(z  +  H)) 

w  =  a - - 

co  cosh  (kH) 

cosh  (k(z  +  H)) 
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cosh(/W) 
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(18) 


It  can  be  seen  from  (18)  that  the  trajectories  of  the  fluid  particles  are  elliptical. 

There  are  several  simplifications  that  may  be  made  to  the  above-derived  expressions  for  the  cases  of 
shallow  (long)  and  deep  (short)  water.  The  shallow  and  deep  water  ranges  correspond  to  H/L  <  71/10  and  H/L  >  K 
respectively,  and  over  these  ranges  approximate  expressions  may  be  substituted  for  the  hyperbolic  functions  that 
have  been  encountered.  Table  1  summarizes  these  results.  Figure  2  depicts  the  water  particle  velocity  between  long 
and  short  waves.  The  key  points  that  should  be  observed  are  that  shallow  water  waves  are  non-dispersive  and  that 
the  classification  of  shallow  water  depends  on  the  ratio  of  water  depth  to  wavelength. 


Table  1.  Summary  of  progressive,  small  amplitude  wave  equations 


Deep  water  (short)  waves, 
H/L>  1/2, 
kH>  7t 

Intermediate  depth  waves, 
1/20  <  H/L  <_l/2, 

7/10  <kH  <n 

Shallow  water  (long)  waves, 
H/L  <  1/20, 
kH  <  7/10 

In  all  cases,  0=fcc-COf 
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co 
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rs 
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o  = - sin  0 
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2k 

Particle  velocity,  u 
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As  an  example,  consider  a  one-meter  high,  monochromatic  wave  with  a  ten-second  period  (0.1  Hz), 
propagating  in  water  six  meters  deep.  The  deep  water  and  shallow  water  wavelengths  are  156  meters  and  76.7 
meters,  respectively.  The  ratios  of  water  depth  to  wavelength  are  0.038,  for  the  deep  water  case,  and  0.078,  for  the 
shallow  water  case.  Referring  to  Table  1,  it  can  be  seen  that  neither  the  deep  water  nor  shallow  water 
simplifications  may  be  used.  In  fact,  the  water  depth  must  be  reduced  to  2.45  meters  before  the  shallow  water 
equations  are  valid.  Conversely,  if  the  water  depth  remains  at  6  meters,  then  the  wave  period  must  become  greater 
than  15.6  seconds  or  the  wave  frequency  less  than  0.064  Hz.  This  indicates  that  only  low  frequency  waves  typically 
qualify  as  shallow  water  (non-dispersive)  waves  which  becomes  an  important  consideration  in  the  proper  modeling 
of  the  wave  induced  disturbances.  This  is  shown  graphically  in  Figure  3. 

As  a  side  note,  it  should  be  pointed  out  that  the  wave  height  has  no  bearing  in  determining  whether  a  wave 
is  classified  as  long  or  short.  The  wave  height  is  significant  in  order  to  determine  the  subsurface  water  particle 
velocities,  and  when  the  wave  breaks.  As  a  rule  of  thumb,  a  wave  will  typically  begin  to  break  when  the  wave 
amplitude  approximately  equals  the  water  depth,  but  this  is  not  by  any  means  exact,  since  other  factors,  including 
wave  steepness  and  bottom  topography,  influence  wave  breaking. 

The  plane  progressive  wave  described  so  far  is  a  single,  discrete  wave  system,  with  a  prescribed 
monochromatic  component  of  frequency  ft)  and  wavenumber  k,  moving  in  the  positive  x-direction.  More  general 
wave  motions,  which  are  no  longer  monochromatic,  can  be  obtained  by  superimposing  plane  waves  of  different  k 
and  co.  The  elevation  of  the  sea  surface  r\(t)  can  thus  be  described  as  the  superposition  of  an  infinite  number  of 
sinusoids  of  the  form: 


11  {t )  =  E 1 an  COS(knX  -  CO,/  +  <|>„  )  =  • 


n= 1 


n=  1 


(19) 


The  equations  for  the  horizontal  and  vertical  velocities,  accelerations  and  the  dynamic  pressure,  respectively,  are: 


n=l  L 


ft),,  cosh  kn  (H  +  z ) 


sinh  k„H 


(20) 


w(t) = £ 


n=l  L 


ft),,  sinh  £,,(//  +  z) 
sinh  k,,H 


hi. 


(21) 


P(t)  =  PS^ 


n=l  L 


ft),,  cosh  kn  (77  +  z) 


cosh  k„H 


h„  -  p  gz 


(22) 


where  the  parameters  in  ( 19)-(22)  are  identical  to  those  in  the  monochromatic  case. 

Wave  Modeling  Using  Recursive  Methods 

A  typical  experiment  in  system  identification  consists  in  recording  a  set  of  input/output  data  and  fitting  a 
parametric  model.  In  discrete  time,  the  first  attempt  is  to  fit  a  Linear  Difference  Equation  (LDE)  model  of  the  form 

y(t)  +  aly(t-l)  +  ...  +  any(t-n)  =  btu(t-])  +  ...  +  bnu(t-n)+£(t)  (23) 


with  l  e  Z  denoting  the  discrete  time  index,  and  E(t)  an  error  term,  accounting  for  the  fact  that  the  data  never 
matches  the  model  exactly.  The  problem  is  to  estimate  the  parameter  vector 

0  =[a1,...,an,b1,...,bnf  (24) 

from  a  set  of  data.  In  this  section  we  will  present  the  general  procedure  to  estimate  the  parameters  associated  with  a 
discrete  linear  model. 

A  way  of  writing  (23)  is  in  regression  form  as 
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y(t)=Ht-lfe  +e(0 


(25) 


with 


4> (r  —  1 )  =  [ —  }’(t  - 1),...,— y{t  -  n),u(t  - 1 u(t  —  n)Y  (26) 

being  a  sliding  window  of  input  and  output  data.  Now  the  problem  is  to  determine  a  technique  to  compute  an 
estimate  of  the  parameter  vector  9  on  the  basis  of  the  data  set  Z  =  {u(0),...,U(N),y(0),...,y(N)},  if  it  is  assumed 
that  N  data  points  are  collected. 

If  the  problem  is  cast  in  a  probabilistic  framework,  a  probability  density  for  the  data  set  Z  given  the  model 
9  can  be  written  as 


P(Z|9 )  =  Pr(e  (0), . . .  ,e  (AO)  (27) 

with  e(t)  =  y(t)—  §T(t  —  7)9  .  If  we  assume  the  disturbance  term  sequence  to  be  gaussian  and  white,  then  the 
density  on  the  right  hand  side  becomes 

N 

Pr(£  (0), . . . ,  £  (AO)  =  J~J  Pr(£  (f))  (28) 

r=0 

with 

£  2 

Pr(e )  =  — - — e  20  2  .  (29) 

27tC7 

As  a  consequence  it  can  be  seen  that 

'  .£!?«>  o' (7  1)0  2 

Pr(Zp)  =  Ce  20  .  (30) 

As  can  be  seen,  minimizing  the  summation  in  the  exponential  can  maximize  this  probability.  So  if  the  estimate  of 
the  parameters  is  defined  as  the  vector  that  minimizes  the  probability,  as 

0  =  arg  mine  Pr(  Z|9 )  (31) 

then  it  can  be  computed  as  a  least  squares  solution 

M  , 

9‘  =  argmine  ^|y(?)-<|)rW9|  .  (32) 

t=\ 


The  solution  can  be  obtained  using  standard  techniques,  by  writing  9  as  the  least  squares  solution  of  the  system  of 
equations 
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or,  compactly 


y  =  $0.  (34) 

The  solution  given  by  the  pseudoinverse  becomes 

0  (35) 


which  can  be  written  as 


0  = 


-lr 


M  “  M  “ 

t= 1  t= 1 


(36) 


This  is  true  provided  the  error  sequence  e( I)  is  white  and  gaussian.  If  the  error  is  not  white,  then  the  estimate  of  the 
parameter  vector  0  is  biased. 

Now  that  we  have  set  the  foundation,  the  problem  that  arises  is  how  do  we  estimate  the  parameters 
associated  with  a  transfer  function  that  will  properly  represent  a  model  of  the  seaway  in  a  recursive  fashion.  If 

0  (f)  is  the  estimate  of  the  parameters  at  time  f,  the  goal  is  to  compute  the  successive  estimate  0  (t  +  I )  by  updating 

0  (f )  using  the  latest  observations.  It  turns  out  that  this  problem  can  be  put  in  a  very  nice  framework  that  makes  use 
of  the  considerations  on  the  Kalman  Filter. 

The  AR  model,  with  numerator  equal  to  one,  can  be  written  as 

y(t)=Q>(t-l)TQ+e(t)  (37) 

with  e(t)  a  white  noise  sequence.  If  it  is  assume  that  0  is  constant,  (37)  can  be  written  in  state  space  form,  where  the 
state  is  the  parameter  vector  itself,  as 


0(f  +  l)=0(O 


This  is  just  a  particular  case  of  the  stochastic  state  space  model  with  A  being  the  identity  matrix,  B=0  and  C=(|>  ( ?  J  7  . 
This  calls  for  a  Kalman  Filter  approach  for  the  estimation  of  0(f),  which  leads  to 


0(f  +  l)=0(f)  +  - 


P(fXKf-l) 


X2  +<Kf-i)7>(f)<Kf-i) 


(y(f)-0(f)'  )4>(f-l) 


P(t  +  l)  =  P(t) 


p(t)<b(t-n\>(t-pT  p(t) 
X2  +<\>(t-l)TP(t)i\)(t-l) 


(39) 


with  X2  =  E{|<?(0|j.  Clearly  this  parameter  X  is  never  known,  and  the  need  to  know  it  can  be  eliminated  by  proper 
normalization.  Dividing  the  numerator  and  denominator  of  (39)  by  X2  the  following  algorithm  is  obtained; 


0(f  +  l)  =0(f)  +  - 


P(tW~  i) 


1+<K*-1)  ^W<K*-1) 


-(y(t)-m'  )4>(f-l) 


P(t+l)=p(t)~ 


p(tm-im-ir  Pit) 
i+ty(t-i)T  P(tMt-i) 


(40) 


which  can  be  easily  verified  by  setting  P(t)  =  P (t)  /  A, '  .  This  algorithm  is  called  the  Recursive  Least  Squares  (RLS) 
algorithm. 

Using  the  above  developed  RLS  algorithm  the  parameters  of  an  Auto  Regressive  model  of  the  form 

w(t)  +  djw(t  — 1)  + ...  +  dNw(t  —  N)  =  e(t) ,  (41) 

with  w(t)  representing  the  sea  surface  elevation  due  to  wave  action,  and  e(t),  a  zero  mean,  white  noise  otherwise 
know  as  the  innovation,  may  be  developed.  The  main  advantage  of  the  AR  model  is  the  fact  that  the  parameter 
estimation  is  a  linear  operation.  In  fact,  if  the  numerator  pertinent  to  the  noise  term  is  not  one,  then  the  error  state 
model  (38)  is  not  white  and  we  cannot  apply  Kalman  Filtering  techniques  directly.  In  the  general  case,  the  problem 
is  nonlinear  in  the  parameters  and  the  Extended  Kalman  Filter  must  be  applied. 

The  problem  that  now  arises  is  to  determine  the  order  (AO  of  the  model  to  properly  reflect  the  actual 
frequency  spectrum  of  the  time  series  w(t)  while  keeping  the  complexity  if  the  model  at  a  minimum.  By  computing 
the  covariance  of  the  innovation  as  the  order  of  the  model  is  increased  it  can  be  shown  that  there  is  a  value  to  which 
the  covariance  converges.  Using  surface  wave  elevation  data  obtained  in  Monterey  Bay,  from  a  Waverider 
measurement  buoy,  an  Auto  Regressive  model  of  various  orders  was  determined  and  the  covariance  of  the 
innovation  compared.  The  results  of  this  analysis  are  shown  is  Figure  4. 

As  can  be  seen,  the  covariance  begins  to  flatten  out  around  an  eighth  order  model.  This  would  have  one 
believe  that  the  correct  order  model  to  choose  would  be  an  eighth  order  model.  However,  as  Figure  5  indicates,  an 
eighth  order  model  fails  to  accurately  reflect  the  actual  spectrum.  Figures  6  through  8  show  that  as  the  order  of  the 
model  is  increased,  the  error  in  the  matching  of  the  actual  spectrum  decreases.  It  is  not  until  a  100th  order  model  that 
the  desired  spectrum  is  actually  realized.  Figure  9. 

There  are  those  that  will  argue  that  the  energy  associated  with  the  low  and  high  frequency  modes  is 
minimal,  and  that  proper  modeling  of  those  modes  is  not  important.  However,  as  was  discussed  earlier,  it  is  the  low 
frequency  modes  of  the  wave  train  that  are  considered  shallow  water,  non-dispersive  waves.  It  is  these  modes  that 
will  have  the  most  effect  on  the  submerged  vehicle  trying  to  maintain  station.  Therefore,  as  in  any  design,  trade-offs 
must  be  made  based  on  analysis,  as  to  how  accurate  the  wave  disturbance  model  must  be  compared  to  the  resulting 
controller  complexity. 

Time  Series  Transformations  Of  Surface  Wave  Records 


The  simplest  method  of  obtaining  information  about  the  wave  disturbances  in  a  particular  operating  area  is 
by  surface  elevation  measurements  from  a  wave  buoy.  Therefore,  it  is  necessary  to  be  able  to  transform  the  wave 
elevation  record  into  a  water  particle  velocity  record  at  the  desired  vehicle  operating  depth  in  order  to  determine  the 
proper  disturbance  model.  There  are  two  approaches  to  this  problem.  The  first  uses  spectral  analysis  while  the 
second  uses  Fourier  analysis. 

Using  the  spectral  analysis  approach,  the  procedure  is  to  first  compute  the  power  spectral  density,  then  use 
the  appropriate  modifier  and  finally  convert  the  spectral  density  back  to  time  domain.  Equation  (42)  shows  this 
relationship. 


5„„(0))  =  |H(co)|2ST1T1((0)  (42) 

Although  the  resulting  time  series  reflects  the  proper  magnitude  of  the  water  particle  velocities,  the  disadvantage  to 
this  method  is  that  all  phase  information  is  lost,  and  the  resulting  subsurface  velocity  time  series  does  not  reflect  the 
motion  caused  by  the  surface  elevation  times  series. 

Through  the  use  of  Fourier  analysis  this  disadvantage  is  overcome.  In  this  method,  a  FFT  of  the  surface 
wave  record  is  taken,  then  the  Fourier  coefficients  are  multiplied  by  the  appropriate  value  of  the  modifier  for  each 
frequency  component.  Once  this  has  been  completed,  then  an  inverse  FFT  is  performed.  The  resulting  times  series 
now  reflects  both  the  proper  amplitude  and  phase  relation. 


Figures  10-12  demonstrate  this  Fourier  analysis  translation  procedure  for  a  wave  elevation  time  series 
recorded  in  30  feet  (9.14  m)  of  water,  and  translated  to  a  subsurface  water  velocity  record  at  a  depth  of  20  feet  (6.1 
m)  in  Monterey  Bay,  CA.  It  will  be  this  translated  velocity  record  that  will  be  used  during  the  control  law  design  and 
vehicle  simulations  later  in  this  paper.  The  maximum  velocity  as  shown  in  Figure  12  is  approximately  1-knot 
(0.5144  m/s).  It  is  interesting  to  note,  see  Figure  13,  that  in  addition  to  the  dominate  peak  frequency  at  0.12  Flz 
, (caused  by  the  wind),  there  are  two  additional  lower  frequency  swell  generated  peaks  at  0.05  Flz  and  0.07  Flz. 
Referring  back  to  Figure  4  and  Table  1,  it  can  be  seen  that  for  this  water  depth  that  only  the  0.05  Flz  frequency 
component  may  be  treated  as  a  shallow  water  wave. 

AUV  Surge  Model  Development 

The  equations  of  motion  (EOM)  based  on  constant  hydrodynamic  coefficients  for  an  underwater  vehicle 
are  well  known  and  documented  in  the  literature  (Flealey,  1993).  In  general,  the  hydrodynamic  force  on  a  submerged 
vehicle  depends  on  the  relative  velocity  and  acceleration  between  the  water  particles  and  the  vehicle.  The  general 
case  of  longitudinal  (surge)  motion  in  the  indirection  written  in  compact  form  is 

(m  +  Xu)ur  +  Xmlur  \ur\  =  Xprop{t)  (43) 

where 

Ur  =  (U~U  fluid)  (44) 

and  Xpmp(t)  represents  the  propulsion  force.  The  propulsion  force  can  also  be  modeled  in  the  form 

-1  (3  I, 

Xprop  -  —Xprop  +~n\n\  (45) 

where  x  and  (3  represent  a  time  constant  and  a  thrust  parameter,  respectively.  Combining  (43)  and  (44)  and  adding 
the  kinematic  relations  for  a  vehicle  constrained  to  longitudinal  motion,  the  following  system  of  equations  is 
developed. 


X  -  ur  +  Ucx  +  Uwx 
ur=aur\ur\  +  Fprop  , 


p  -7  (3  I. 

FproP=  —  F+—n\n\ 
11  x  x 


(45) 


where  because  relative  velocities  are  used  in  (45),  the  steady  state  current  and  fluctuating  velocity  caused  by  waves 
can  be  added  directly  to  the  kinematic  equation.  It  should  be  pointed  out,  that  in  the  form  used  in  (45),  F  becomes  a 
generalized  force  with  unit  the  same  as  an  acceleration  vice  a  direct  thrust  value. 

With  the  development  of  (45)  complete,  it  is  desirable  to  linearize  this  system  of  equations  about  some 
nominal  operating  point  to  simplify  the  controller  design.  This  nominal  operating  point  can  be  taken  to  be  a  variety 
of  conditions  depending  on  the  analysis  and  control  development  to  be  performed.  In  our  case,  since  we  are 
interested  with  station  keeping,  we  chose  to  linearize  about  the  steady  state  solution  to  (45)  in  the  presence  of  a 
steady  current  and  obtain 


=  -u, 


F  —  (i  1 !  z 
r  prop, o  'JjlJ  c. 


aU: 


(46) 


as  the  nominal  operating  point.  Using  a  standard  Taylor  series  linearization,  the  linearized  system  of  equations,  in 
state  space  from,  can  then  be  written  as 


9 


X 

0 

1 

0 

~x~ 

0 

~r 

«r 

= 

0 

2a  sgn  (u0)uo 

1 

ur 

+ 

0 

2—sgn(n0)n0 

n  + 

0 

F 

0 

0 

-1 

F 

0 

T 

T 

"l 

0 

0~ 

'X' 

y  = 

0 

1 

0 

ur 

0 

0 

1 

F 

(47) 


assuming  that  all  states  are  measurable.  This  may  be  a  little  pretentious,  but  it  is  possible  to  obtain  X  and  ur  from  a 
doppler  sonar  with  bottom  lock,  and  it  would  be  possible  to  get  thmst  measurements  directly  from  the  propulsion 
system  if  it  had  strain  gauges  embedded  in  the  propeller. 

If  it  is  assumed  that  the  vehicle  will  be  operating  in  a  1-knot  (0.5144  m)  steady  current,  and  the  parameters 
a,  p.  and  x  are  available,  then  equation  (47)  can  be  evaluated  numerically.  As  Marco  (Marco,  1998)  has  pointed  out, 
through  system  identification,  it  is  possible  to  obtain  these  values.  The  values  obtained  for  PHOENIX  and  reported 
at  IARP  '98  are 


a  =  -4.036  m  1 

P  =3.962  xl0~3  m/rev2  .  (48) 

x  =3.14  s 

Since  the  time  of  Marco's  experiments,  the  PHOENIX  has  had  its  propulsion  system  upgraded,  to  include 
larger  brushless  DC  motors  for  propulsion,  increased  diameter/pitch  propellers  and  ducted  shrouds.  From  initial 
at-sea  experiments,  the  vehicle,  due  to  these  system  upgrades,  is  now  capable  of  approximately  4-knots  (~2  m/s)  at 
800  revolutions  per  minute.  Based  on  these  upgrades,  a  more  realistic  set  of  values  is 

a  =  -4.036  m_; 

P  =1.034  m/rev2  ,  (49) 

T  =1.15  s 

which  will  be  used  throughout  the  remaining  design  and  simulation. 

Since  we  will  be  using  parameters  obtained  from  a  discrete  filter,  and  desire  to  implement  the  designed 
control  law  into  PHOENIX,  the  state  space  equations  must  be  converted  into  a  discrete  form.  Using  a  standard  Euler 
discretization  (47)  can  be  represented  in  the  discrete  form  as 


X(k  +  1 )" 

0 

dt 

0 

'*(*)' 

0 

dt 

ur(k  + 1 ) 

= 

0  1 

-  2asgn{u0)uodt 

dt 

ur(k ) 

+ 

0 

2—sgn(n0)n0dt 

T 

n(k)  + 

0 

F(k  +  1) 

0 

0  1- 

-1  , 

- dt 

T 

F{k)_ 

0 

y(k)  = 


1  0  0 
0  10 
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X(k) 

u,.(k) 

F(k) 


Uwx(k) 


(50) 
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Optimal  Predictive  Controller  Development 

Using  standard  optimal  control  techniques  the  solution  for  the  optimal  (LQR)  controller  can  be  found  as 

n  =  -R~‘  BT  Sx  (51) 

where  S  is  found  by  solving  the  steady  state  algebraic  Riccati  equation  (ARE)  for  the  positive  definite  S, 

AtS  +  SA-SBR~1BtS+Q=0.  (52) 

Here  Q  is  the  weighting  matrix  on  the  state  error,  and  R  is  the  scalar,  since  this  is  a  single  input  system,  that  envokes 
a  penalty  against  the  control  effort.  The  LQR  approach  will  always  yield  a  stable  system,  as  long  as  the  Riccati 
equation  provides  a  positive  definite  solution  matrix  S. 

Figure  14  shows  graphically  the  relation  between  the  level  of  control  input  to  the  level  of  disturbance 
rejection  for  the  standard  optimal  control  solution,  without  a  disturbance  model  included  in  the  control  law  design. 
This  analysis  can  give  a  "feel"  for  how  tight  a  control  law  must  be  provided  to  achieve  a  reasonable  disturbance 
rejection.  Using  the  control  law  design  that  provides  the  greatest  disturbance  rejection,  it  can  be  seen  in  Figures  15 
and  16  that  a  reasonable  position  is  maintained  within  the  propulsion  system  limitations.  To  determine  how  well  the 
propulsion  system  can  respond  in  the  presence  of  wave  induced  disturbances,  the  velocity  time  series  was  scaled  to 

a  maximum  value  of  1.5  m/s.  At  this  level  of  disturbance,  the  error  in  maintaining  position  is  +/-  0.4  meters.  The 

required  control  input,  to  maintain  this  position,  is  greater  than  the  input  available  from  the  propulsion  system. 
Figures  17  and  18.  If  the  control  law  is  modified  to  limit  the  controller  to  provide  an  input  within  the  propulsion 
system  design,  then  the  optimal  position  control,  for  this  disturbance,  is  shown  in  Figures  19  and  20.  The  problem 
that  now  must  be  solved  is  how  to  achieve  better  performance.  We  propose,  that  by  embedding  a  predictor  of  the 
disturbance,  into  the  control  system  design,  we  can  achieve  improved  station  keeping. 

The  AR  model  of  the  wave  disturbance  may  be  written  in  state  space  form  as 

Xw{k+l)  =  AwXw{k)  +  Bwv{k) 
yw(k)  =  CwXw(k ) 

with  the  only  measurement  that  is  available  is  the  wave  velocity  at  current  time  k  obtained  from  the  Doppler.  If  the 
system  is  forward  time  shifted,  then  information  about  the  future  disturbance  is  available  based  on  current 
measurements.  Equation  (54)  represents  this; 

wft  +  N)  +  djw(t  +  N  - 1)  + . . .  +  dNwft)  =  eft) .  (54) 

This  is  possible  since  the  basis  of  the  model  is  that  eft)  is  white  noise.  Now  by  augmenting  the  vehicle  state 
equations  with  the  disturbance  state  equations,  a  new  control  law  may  be  developed  using  the  predicted  disturbance 
states.  To  do  this  a  new  state  vector  is  defined  as 


X 


aug 


X 


vehicle 


(55) 


where  the  disturbance  states  are  given  as 


Xw=[Xw(k  +  N-l),...,Xw(k  +  l)]T  ,  (56) 

then  using  the  separation  principle,  the  control  law  is  designed  assuming  all  states  are  measurable.  Therefore,  as  in 
the  previous  optimal  control  discussion,  the  ARE  is  solved  to  obtain  the  appropriate  gains  for  the  augmented  system. 
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This  augmented  system  is 


Xaug(k  +  1)  = 


A  FC 

^vehicle  1  w 


xaug(k)+ 


B 


vehicle 


n(k )  + 


0 

B„, 


v  (k) . 


(57) 


With  the  control  law  designed  the  estimator/predictor  must  be  designed.  Using  optimal  estimation  theory,  a 
estimator/predictor  of  the  form 


XJk  +  T)  =  (Aw  -  LCJXw(k )  +  Lw(k ) ,  (58) 

where  w(k)  is  the  current  disturbance  measurement,  was  developed. 

To  display  how  well  this  procedure  can  work  if  an  accurate  model  of  the  disturbance  is  available,  consider 
the  simple  case  of  a  sine  wave  with  only  one  frequency  component  as  the  wave  disturbance.  The  precise  model  of 
this  disturbance  is  known.  When  this  model  is  embedded  in  the  control  system  design,  perfect  cancellation  of  the 
wave  disturbance  effects  on  station  keeping  may  be  obtained.  These  results  are  displayed  in  Figures  21  and  22.  Now 
these  results  are  for  demonstration  purposes  only,  and  we  do  not  expect  to  get  perfect  cancellation  of  the  wave 
disturbances  since  an  exact  model  of  the  random  sea  is  not  available. 

Using  the  same  weighting  values  that  went  into  the  design  of  the  control  law  that  was  used  in  the 
simulation  of  the  results  presented  in  Figures  19  and  20,  improved  results  are  achieved  from  this  method  when  an 
eighth  order  AR  model  is  used  to  form  the  predictor,  see  Figures  23  and  24.  As  can  be  seen,  there  is  a  400% 
improvement  in  station  keeping,  and  the  control  input  requirements  are  significantly  less.  This  reduction  in  control 
effort  is  particularly  important  given  the  fact  that  power  consumption  is  the  downfall  of  AUVs. 

Conclusions 


This  paper  has  shown  that  an  Auto  Regressive  model  provides  a  convenient  means  of  matching  the  wave 
induced  disturbance  dynamics.  It  has  been  concluded  that  a  digital  estimator/predictor  derived  from  a  low  order  AR 
model  provides  adequate  prediction  of  the  disturbances  acting  on  the  submerged  vehicle.  What  cannot  be  concluded 
from  this  particular  work  is  the  degree  of  optimality  to  which  the  embedded  predictor  provides,  since  this  analysis 
was  done  under  ideal  conditions,  i.e.  no  noise  and  all  states  available.  Thus  far,  based  on  previous  studies  in  this 
area,  the  work  conducted  herein  has  indicated  a  promising  direction  for  improved  prediction  of  wave  disturbances 
on  submerged  vehicles  and  the  resulting  control  required  to  allow  the  vehicle  to  maintain  station. 

The  estimator  developed  is  merely  a  baseline  version  for  predicting  the  disturbance.  That  is,  since  the  AR 
model  from  which  the  estimator  was  developed  was  assumed  to  have  linear  time  invariant  parameters,  the  estimator 
may  not  accurately  reflect  the  actual  dynamics  of  the  disturbance. 

Lastly,  it  should  be  noted  that  the  optimal  predictive  control  law  developed  provides  excellent  disturbance 
rejection  in  the  presence  of  wave  induced  disturbances. 

Recommendations 


First,  it  is  recommended  that  further  investigations  be  made  into  the  modeling  of  wave  induced 
disturbances.  The  ability  to  identify  the  predominate  frequency  components  and  use  a  frequency  weighted  non¬ 
linear  RLS  algorithm  to  determine  the  disturbance  model  parameters  used  in  the  controller  design  will  improve  the 
ability  of  the  vehicle  to  maneuver  and  hold  station  in  the  presence  of  wave  induced  disturbances. 

Second,  the  parameters  of  the  disturbance  model  should  be  considered  to  be  time  varying,  and  an  adaptive 
algorithm  developed  to  increase  the  accuracy  of  the  disturbance  model  used  in  the  control  law  design.  This 
parameter  identification  can  be  done  as  a  background  process,  with  the  resulting  parameters  passed  to  the  controller. 

Lastly,  models  of  vehicle  sensor  performance  in  shallow  water  must  be  developed  to  allow  for  a  detailed 
analysis  to  determine  how  and  which  sensor(s)  should  be  used  to  measure  the  wave  induced  disturbances,  and  how 
accurate  these  measurements  will  be. 
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Figure  2.  Comparison  Of  Water  Particle  Velocities  For  Different  Types  Of  Waves 


I 


Normalized  Power  Spectral  Density 


0.9 


1 


Figure  6.  16th  Order  AR  Model 
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Figure  9. 100th  Order  AR  Model 
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Figure  10.  Surface  Elevation  Time  Series 
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Figure  11.  Modifier  For  Translation  Of  A  Surface  Elevation  Time  Series  To  A  Subsurface  Velocity  Time 
Series  For  A  Water  Depth  Of  9.14  Meters  (30  Feet)  And  An  Operating  Depth  Of  6.1  Meters  (20  Feet) 


Figure  12.  Subsurface  Velocity  Time  Series  At  Operating  Depth  Of  6.1  Meters  (20  Feet) 
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Figure  13.  Power  Spectral  Density  Of  Subsurface  Velocity  Record 


Figure  14.  Comparison  Of  Control  Input  Covariance  To  Normalized  Vehicle  Position  Covariance 
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Figure  15.  Results  Of  Vehicle  Station  Keeping  Using  Optimal  Control  Solution  In  The  Presence  Of  A 

1-Kt  (0.5144  ni/s)  Wave  Induced  Disturbance 


Figure  16.  Propeller  RPMs  From  Optimal  Control  Solution  In  The  Presence  Of  A  1-Kt  (0.5144  m/s)  Wave 

Induced  Disturbance 


21 


Time  (sec) 


Time  (sec) 


Figure  18.  Propeller  RPMs  From  Optimal  Control  Solution  In  The  Presence  Of  A  3-Kt  (1.5  m/s)  Wave 

Induced  Disturbance 
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Figure  19.  Results  Of  Vehicle  Station  Keeping  Using  Optimal  Control  Solution  In  The  Presence  Of  A 
3-Kt  (1.5  m/s)  Wave  Induced  Disturbance  With  A  Control  Input  Within  Operating  Limits 


Figure  20.  Propeller  RPMs  From  Optimal  Control  Solution  In  The  Presence  Of  A  3-Kt  (1.5  m/s)  Wave 
Induced  Disturbance  With  A  Control  Input  Within  Operating  Limits 
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Figure  21.  Results  Of  Vehicle  Station  Keeping  Using  Predictive  Optimal  Control  Solution  In  The  Presence 

Of  A  2-Kt  (1  m/s)  Sine  Wave  Disturbance 
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Figure  22.  Propeller  RPMs  From  Optimal  Predictive  Control  Solution  In  The  Presence  Of  A  2-Kt  (1  m/s) 

Sine  Wave  Disturbance 


Vehicle  Position  (meters)  Amplitude 


0  5  10  15  20  25  30 


Time 

Figure  23.  Vehicle  Relative  Velocity  And  Wave  Induced  Velocity  vs.  Time 
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Figure  24.  Results  Of  Vehicle  Station  Keeping  Using  Embedded  Disturbance  Model  And  Identical  Weighting 
Values  As  Those  Associated  With  The  Results  In  Figure  19 


